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Abstract 

A two-dimensional nonlinear Schrodinger lattice with nonlinear coupling, modelling 
a square array of weakly coupled linear optical waveguides embedded in a nonlin- 
ear Kerr material, is studied. We find that despite a vanishing energy difference 
(Peierls-Nabarro barrier) of fundamental stationary modes the mobility of local- 
ized excitations is very poor. This is attributed to a large separation in parameter 
space of the bifurcation points of the involved stationary modes. At these points 
the stability of the fundamental modes is changed and an asymmetric intermediate 
solution appears that connects the points. The control of the power flow across the 
array when excited with plane waves is also addressed and shown to exhibit great 
flexibility that may lead to applications for power-coupling devices. In certain pa- 
rameter regimes, the direction of a stable propagating plane-wave current is shown 
to be continuously tunable by amplitude variation (with fixed phase gradient) . More 
exotic effects of the nonlinear coupling terms like compact discrete breathers and 
vortices, and stationary complex modes with non-trivial phase relations are also 
briefly discussed. Regimes of dynamical linear stability are found for all these types 
of solutions. 
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1 Introduction 



The use, in optics, of materials with nonhnear responses to external fields 
has over the past years provided a wealth of applications and prospects for 
all-optical communication, where light manipulates light itself [1]. Especially, 
nonlinear photonic structures endowed with a spatially periodic refractive in- 
dex provide new means to control the flow of light in the form of localized 
wave packets, not possible in bulk media [2] . The experimental observation of 
discrete localization in one and two dimensions, in arrays of coupled waveg- 
uides experiencing a cubic (Kerr) nonlinearity [3,4,5,6] or a saturable non- 
linearity through the photovoltaic effect in a photorefractive medium [7] as 
well as in optically induced nonlinear photonic lattices [8,9,10,11], shows sig- 
nificant progress in this direction. The theoretical foundation and predictions 
of these phenomena can, to good approximation, be made by models of the 
discrete nonlinear Schrodinger (DNLS) type [6,12,13,14,15,16] derived within 
coupled-mode theory [17] or the tight-binding approximation using Wannier 
function expansion [18]. The latter technique is also applicable for models of 
Bose-Einstein condensates trapped in deep optical lattices. In these contexts 
the issue of mobility of localized packets of energy is highly interesting as it 
provides insights as how to control optical signals, e.g., for applications in 
multiport switching. 

As is well known [19], travelling localized solutions, or solitons, quite generally 
exist in continuum systems as a result of the competition between nonlinear- 
ity and dispersion. When the continuous symmetry is broken, as is the case 
when considering a nonlinear discrete system, the generic solution is instead 
a non-moving time-periodic one, called a discrete breather (DB) or intrinsic 
localized mode (see, e.g., [20] for a review). The presence of a lattice leads 
to pinning of localized solutions, as this underlying structure effectively intro- 
duces a periodic potential for travelling excitations due to the difference in 
energy between stationary solutions localized on a lattice site and solutions 
localized between lattice sites [21,22]. This energy difference is referred to as a 
Peierls-Nabarro (PN) energy barrier and is the reason that narrow travelling 
excitations get trapped in a lattice. Mobility of a highly locahzed stationary 
mode can be induced if an appropriate perturbation is added to overcome the 
potential barrier, but usually some energy is lost to radiation and the result 
is trapping [23,24]. However, excitations that are broad, in comparison with 
the lattice spacing, can still propagate as the effects of the PN-barrier in this 
case is small and the dynamics is well approximated by the continuum limit 
of the discrete equation. This is provided when the continuum equation has 
stable localized travelling solutions, which, e.g., is the case for the ID cubic 
nonlinear Schrodinger equation but not for the 2D (cubic) equation [25]. 

Notable exceptions from the pinning behaviour are integrable lattice equa- 
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tions, i.e., equations possessing an infinite number of conserved quantities, 
like tlie Ablowitz-Ladik (AL) model which is a discretization of the nonlinear 
Schrodinger equation retaining this property [26]. These models have analytic 
exponentially localized travelling solutions and completely lack PN-barrier, 
which accounts for the excellent mobility. Unfortunately, in contrast to their 
many interesting mathematical properties, they seem to have very limited di- 
rect applications for real systems. Thus other routes to mobility in discrete 
systems are sought, although these models are conveniently used as the start- 
ing point of perturbative approaches. In [27] it is numerically shown how a 
mobile excitation develops a resonant tail when the governing equation devi- 
ates from the ID integrable AL-model. Moreover, the amplitude of the tail 
grows when the equation is further away from the integrable limit, or alter- 
natively, as the PN-barrier increases from zero. It is inferred that the effect of 
pinning is overcome by a periodic exchange of energy between the localized 
core and the oscillating tail as the excitation moves through the lattice and the 
periodic PN-potential. This result can also be extended to a 2D setting [28]. 
However, there is a limit as to how far the mobile solutions can be continued 
away from the integrable limit, with broad solutions being more persistent as 
can be expected [27]. Hence, we can expect that mobility of localized modes 
is possible when the PN-barrier is small, and the abihty to tailor the size of 
the PN-barrier may allow for a control of the mobility of narrow solutions in 
a lattice. 

Recent research has focused on ways to minimize the PN-barrier, which in 
optics applications means moving beyond models based on the cubic DNLS 
equation where the PN-barrier is an increasing function of the power of the 
stationary modes [4]. Similar to the problem of moving kinks in Klein-Gordon 
lattices [29] the introduction of competing nonlincaritics has proven fruitful. In 
DNLS equations with saturable nonlinearity the vanishing of the PN-barrier 
at certain points and a subsequent enhanced mobility has been demonstrated 
in both ID [30] and 2D [31], and the same apply for a ID DNLS model with 
nonlinear interactions [32]. As these models are derived from real physical 
settings they are interesting from an application point of view. Note also that 
there is a large class of DNLS models [33] with an absence of PN-barrier (not 
only zero at isolated points), different from the more physical models, as the 
former are special discretizations of the nonlinear Schrodinger model that in 
some sense preserve a translational invariant [34,35]. Amongst them is the 
AL-lattice and other integrable equations. 

A problem related to the vanishing of the PN-barrier at specific points is that 
of stability inversion between stationary solutions centred on a lattice point 

and between lattice points, respectively, which can be understood from the 
energy minimum principle. The solution with lowest energy at given power 
will be stable, and the change of sign of the energy difference is also related to 
an exchange of stability of the lowest energy modes. This can be observed for 
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the above models [30,31,32,36,37], as well as for the DNLS model with com- 
peting cubic-quintic nonlinearity [38] and more general lattice models [39,40]. 
In connection with the stability inversion there also exists an intermediate 
asymmetric solution interpolating between the stationary solutions centred 
on and in-between sites and bifurcating with these at the points where they 
change their stability properties [31,32,38,39,40]. For the ID saturable DNLS 
model this solution can even be found analytically near the first zero of the 
PN-barrier [41]. From this viewpoint, the mobility of a highly localized mode 
may be thought of as a transformation between stationary on- and off-site so- 
lutions, via the intermediate solution. Hence a more accurate definition of the 
PN-barrier would need to include also the intermediate solution as this may 
be the mode with the highest energy. Therefore, we will distinguish between 
the energy difference between two fundamental modes and the PN-barrier re- 
lating to the energy of all involved stationary modes. Also, in 2D, one would 
in principle need to consider modes corresponding to different directions of 
possible propagation [31]. 

Though localized modes and their mobility are interesting for beam steering 
in waveguide arrays, transport of energy in the transversal direction (from 
waveguide to waveguide) can be achieved also by a collective excitation of the 
waveguides. With a plane wave propagating in the array, the same power will 
be transported along all waveguides, but there will also be a transfer of power 
between the waveguides depending on the wave number (phase gradient) of 
the plane wave. In [42] it is demonstrated that the direction and magnitude of 
the transversal flow of power in a ID array with nonlinear coupling are quite 
versatile under the variation of the amplitude and wave number of the plane 
waves. Here the discussion is extended to 2D and will show an even greater 
flexibility. 

The outline of the Paper is as follows: In Sec. 2 we present a model for opti- 
cal waveguides coupled in two spatial dimensions and embedded in a nonlin- 
ear Kerr material, discuss its relevance and introduce some of the important 
properties of the model. We should stress that our model does not attempt to 
describe conflgurations where the waveguides themselves are nonlinear, such 
as those presently experimentally realized in, e.g., [3,4,5,6,7,8,9,10,11], but 
rather correspond to the original proposal of Jensen [12] with the nonlinearity 
residing solely in the embedding material. The stability of the fundamental 
localized modes are investigated in Sec. 3 and their relation to mobility of 
localized excitations is discussed in Sec. 4. In Sec. 5 some special solutions, 
like compact and complex modes, arising from a balance of the coupling terms 
in the model are presented. Finally, in Sec. 6 the transport properties of plane 
waves are analyzed before we conclude in Sec. 7. We also include in the end 
of Sec. 7 a brief discussion about some recent related works, which appeared 
in press after the original submission of this manuscript. 
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2 Model 



The model we consider is a direct extension to two spatial dimensions of a 
one-dimensional model for the propagation of the electric field amplitudes in 
coupled optical waveguides embedded in a nonlinear material [32]. The waveg- 
uides arc identical and regularly spaced in the two spatial dimensions, thus 
forming a square lattice. Coupling between the sites of the lattice is mediated 
by an evanescent field overlap of the modes (single-mode propagation) in the 
respective waveguides and it is assumed that the modes decay sufficiently fast 
outside the waveguides to motivate only nearest-neighbour coupling. Further, 
the waveguides themselves are assumed to be constructed from a linear mate- 
rial and surrounded by a nonlinear Kerr material. This will tend to strengthen 
the effects of nonlinear coupling with respect to on-site nonlinear effects when 
compared to a system where also the waveguides are nonlinear. This moti- 
vates an extension of the generally employed cubic on-site DNLS model to an 
equation that in two dimensions will read 



i ' = Ql'^n,m + Q2^^n,m + 2(^3^ 

n,m I ^ n,m \ 



dz 



+ 2Q^ 
+ 2Q5 



2\E' Afl\E' -I- \['* Af^f^ 1 

^^n,m^U^n,m| J ^ ^ n,m'-^\^ n,mJ 
2 



(1) 
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where 'ifn,m is the complex amplitude of the electric field, z measures the 
propagation along the waveguide array and the coupling parameters Q1-Q5 
are given by overlap integrals of the modes. The operator A is defined by 

+ '^n,m+i- (Note that the operator A — 4 
is the standard 2D discrete Laplacian.) 

The model (1) is also relevant in other contexts, e.g., as a nonlinear tight- 
binding approximation, with Q4 argued to be negligible, for the time-evolution 
of Bosc-Einstcin condensates (BEC) in a deep periodic optical potential [43,44,45], 
or for = Q4 = as a. rotating- wave approximation to a Fermi-Pasta- 

Ulam chain [21]. In both these one-dimensional models the extension to two 
spatial dimensions is straightforward. One may also in a straightforward way 
consider extensions to other lattice geometries, such as e.g. triangular/hexagonal 
patterns (which in some sense could be considered as the most natural packing 
of waveguides in two dimensions). In this work, we chose to consider only the 
square lattice which is the simplest 2D structure to analyze from a modelling 
point of view, and also for comparison with previously published work in the 
area which mainly dealt with square-lattice configurations. We expect that 
most of the qualitative results obtained here will persist also for different lat- 
tice geometries, although a more detailed discussion on these issues is beyond 
the scope of the present work. 
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2.1 Parameter estimates 



For general waveguide arrays Q2 > and Qj, j = 3,4,5, have the same 
sign as the Kerr index of the nonhnear material. Estimates of the relative 
strengths of the parameters show that the nonlinear coupling terms, Q4 and 
Q5, can be up to the same order of magnitude as the on-site nonlinearity 
for an array of waveguides constructed from AlGaAs with sizes in the tenth 
of micrometer regime and operated with a laser in the infrared (A ~ 1.5/im). 
Details of the derivation of the ID model can be found in [32,46] and a cal- 
culation of the coupling parameters for a slab waveguide array can be found 
in Appendix A of [46] (see in particular Fig. A2 in [46] for an explicit illus- 
tration of realistic values for relative coupling strengths) . When extending the 
model to two dimensions there is also the possibility of coupling in the diag- 
onal direction between next-nearest neighbours. However, for experimentally 
relevant sizes, as given above, of an array of square waveguides the linear di- 
agonal coupling can be estimated to be an order of magnitude smaller than 
the linear direct coupling (see Appendix A2 in [46]). Also, when operated in 
the power regime where the most interesting phenomena occur, the nonlinear 
coupling and on-site terms will be of the same order as the linear coupling, 
i.e. \Q2\ ~ \Qj\ sup{|\E'„ mP}, j = 3,4,5. This can be shown to correspond to 
powers of a few kW [46], well within what is experimentally available. How- 
ever, although these estimates of the parameter values are of relevance when 
considering experimental verifications and applications of the model they will 
not be taken as a restriction in the present paper, since we are also interested 
in the general effects of nonlinear coupling. 

Further, the parameter Qi can be made to vanish by the transformation 
^n,m(^) ^ ^n,m(^)e~''^^^ and has no significant physical role. Note also that 
the staggering transformation ^n,m{z) ^ (— l)""'""^^„,^(z) is equivalent to 
changing the sign of the parameters Q2 and Q5, thus reducing the part of 
parameter space that needs to be investigated to get a complete picture of the 
equation. Further, by rescaling the amplitudes and the time-like variable, two 
additional parameters can be fixed. Here we will fix (^2 = 0.2 and Q3 — 0.5, 
while restricting (^4 > without loss of generality. 



2.2 Conserved quantities 

From Noether's theorem we know that any infinitesimal transformation leav- 
ing the action integral for an equation invariant leads to a quantity that is 
conserved under the evolution governed by that equation [47] . Essentially this 
means that any continuous symmetry of Eq. (1) corresponds to a conserved 
quantity. In particular, the invariance under translations in z will lead to con- 
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servation of the Hamiltonian Ji — Y,n,m'Hn,m, with 



n.ml 



~l~ |^n,m+l| ) 

+ lT/2 ('\T/*2 +\I>*2 )] (2) 

' n.rn\ ri+l,m ' n.m+l/ \ / 

+ 2Q,kJn,m'^r^+l,m{K'!m + ^^'ij 



+ C.C. 



The evolution equation (1), and its complex conjugate, can be obtained from 
the Hamilton equations of motion using the canonical variables ^'n.m and 
i^* , 



(3a) 



"'^ n,m 
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dz 


n,m 


n,m 
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dz 


^ ^ n.m 



-i— ^ = (3b) 



Further, the invariance of Eq. (1) under global phase rotations, i.e., ^n,m ^ 
*n,me'", CK e R, corresponds to the conservation of the norm 

N=Y.Nn,m = Y.\'^n,m?. (4) 

n,m n,m 

Although the Hamiltonian Ti. for all practical purposes constitutes an energy 
functional it has no direct connection to any physical energy for the array 
of waveguides. But from a mathematical viewpoint it can be treated as the 
energy of the system. In particular, for the cubic DNLS model it holds that a 
bounded extremum (a maximum in the present formulation since (^2, Qs > 
and hence Ti. should be regarded as negative energy) of Ti. for given norm 
A/" is Lyapunov stable and may, with right, be called the ground state of the 
system [48]. This ground state is conjectured to be unique up to a global 
phase factor, which is consistent with the absence of stability inversion and 
the non-vanishing PN-barrier in the cubic on-site DNLS model. The norm 
A/", however, has a direct physical interpretation as the total power carried 
along the waveguides, provided that the modes are appropriately normalized. 
In the context of BEC, the conservation of norm corresponds to boson number 
conservation. 

The conservation of Hamiltonian and norm may also conveniently be expressed 
in terms of discrete continuity equations. For this purpose it is instructive to 



make a change to action-angle variables ^'n,m = y-^n,m^ ^Sn,m^ where Nn,m 
and 9n,m are real canonical variables. The corresponding Hamilton equations 
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of motions are 





m 


dz 


dNn,m 




m 


dz 





(5a) 
(5b) 

Now, with T-Ln,m from Eq. (2) and by introducing the phase difference between 
neighbouring sites in the two lattice directions 

<Pn,m — dn,m ~ ^n,m-i Eq. (5b) for the uorm density can be written as 

' '~^n,m »^n— l,m '-^n,m »^n,m— 1 ^" v / 

The norm current density in the n direction is given by 



Jn% = - ^ = '^J J^n,mJ^n+l,m Sin ^n+l,r. 



X 



Q2 + ^QA\lMn,mKi+l,m COS (^„+i . 



(7) 



A similar expression holds for J^i^ = —&Hn,m/ d<pn,m+i in the m direction 
where the roles of n and m are interchanged. Since the norm density Hn,m is 
the power on site (n, m), the continuity equation expresses the obvious fact 
that a change of the power is related to a power transfer to neighbouring 
sites. Hence the physical interpretation of J"^'^^ is the flow of power from site 
(n, m) to site (n + l,m). An equation on the form (6) can also be derived 
for the Hamiltonian density 'Hn,m, with a Hamiltonian flux density X^"2i 
T^^. For the present purposes it is sufficient to know that for a stationary 
monochromatic solution, i.e., Mn,m{z) = Afn,m and 9n,m{z) = 6n,m + Az, the 
relation to the norm current density is T^^}^ = Aj7^'^\. Our discussions on the 
flow of energy in the system (Sec. 6) will thus be made using only the norm 
current density. 



3 Stationary solutions 



Initially we will turn our attention towards localized solutions or, more specifl- 
cally, to time-periodic locahzed modes also known as discrete breathers (DBs) 
or intrinsic localized modes. The existence of DBs has been rigorously proven 
for Hamiltonian lattice systems of arbitrary dimension subject to conditions 
of anharmonicity and non-resonance with linear phonons [49], and extended 
to more general systems in [50] . In the context of coupled waveguides the vari- 
able z will play the role of time. Due to the global phase invariance of the 
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DNLS-type equations the time-periodicity leads to a significant simplification 
of the problem of finding localized solutions. For monochromatic solutions a 
simple gauge transformation to a rotating frame of reference will render these 
solutions stationary, i.e., taking ^n,m{z) = 4>n,m^~^''^'^'^^''^ reduces Eq. (1) to 



n,m I 

+ [2^n,m^(\^l^n,m\') + rn,MlJ 

+ 2Q5 [2|Vn,m|'AV'n,m + ^lm^rn,m 
+ A{tpn,m\'ipn,m\^) =0. 



(8) 



The idea of the existence proofs, and also the foundation for efficient numerical 
techniques to calculate DBs [51,52], is continuation from the anti-continuous 
Unfit, Qj = 0, j = 2,4,5. Solutions are trivial in this uncoupled limit as the 
amplitude of the oscillators at each site can be chosen independently from the 
set ipn,m e {0, ±yj A/2Qz}, where we have made a restriction to real single- 
frequency solutions. From the configuration chosen in the anti-continuous limit 
it is possible to make a classification of the most fundamental localized solu- 
tions of Eq. (8). The set of solutions obtained from continuation to non-zero 
coupling of the configuration with one site excited and all other at rest will 
be called 1-site DBs. Similarly, if two neighbouring sites are excited at the 
ant i- continuous limit we will get symmetric or anti-symmetric 2-site DBs, de- 
pending on the relative phases of the sites. 

Since we ultimately are interested in the mobility of localized modes, we would 
like to compare the Hamiltonian of different stationary solutions. As has been 
described in Sec. 1 (especially for the corresponding ID model [32]), the en- 
hanced mobility of localized modes is connected to the vanishing of a PN- 
barrier measuring the difference in Hamiltonian between fundamental sta- 
tionary modes. The idea is that for a small or vanishing barrier, mobility of a 
highly localized mode may be thought of as a continuous transformation be- 
tween, e.g., stationary 1-sitc and 2-site solutions, put into motion by adding 
energy to overcome the barrier (or pinning). Therefore, numerical continuation 
to non-zero coupling will be performed under the constraint of constant norm. 
As the parameters Q2 and Q3 are fixed from rescalings, the parameters Q4 and 
Qs are varied, leaving the frequency A as a free parameter to be determined 
by the stationary solution. The choice to keep A/" fixed is motivated by the 
fact that it is a dynamical parameter that is conserved under the evolution 
of Eq. (1). Hence, when considering a range of coupling parameters (as will 
shortly be done) , this will facilitate a meaningful comparison of the Hamilto- 
nians over that entire range. This would not be the case had the continuation 
been performed for constant frequency. However, for a fixed frequency and 
varying norm the relevant energy quantity would instead be the grand canon- 
ical free-energy Q = 7i — AJ\f, with A acting as a chemical potential, as used 
in 
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Fig. 1. Region of stability (shaded) for (a) 1-site, (b) symmetric 2-site and (c) 
anti-symmetric 2-site solution for Q2 = 0.2, Qs = 0.5 and A/" = 5; (d) region of exis- 
tence (shaded) of intermediate solutions (unstable) between the stability boundaries 
of the 1-site solution and the respective 2-site solution [symmetric (anti-symmetric) 
in the right (left) region]. Additional small regions of stability, not shown here, can 
be found for the 2-site solutions for relatively small (54- In the cases of multiple solu- 
tions, as for the 1-site case, stability is indicated if at least one solution is stable. See 
also Figs. 3 and 4 for some details. The solid line in (a), (b) and (d) indicates where 
the energy difference, A7^ = Hi — 'H2s-, vanishes for the 1-site and symmetric 2-site 
solutions. The dashed line in (a), (c) and (d) shows the same thing for the 1-site 
and anti-symmetric 2-site solutions, A7Y = 7ii — li.2a- The general phenomenology 
is similar also for other values of the norm M. The size of the lattice is 21 x 21 with 
periodic boundary conditions. 

3.1 Stability 



The linear stability of stationary modes is most simply investigated by study- 
ing the evolution of a small perturbation in the rotating frame of reference. 
Hence, for a given solution we make the ansatz ^n,m{z) — [■0n,m + {^n,m + 
i'rin,m)^^^]^'~^^^^^^^^ in Eq. (1) and arrive at a linear matrix eigenvalue prob- 
lem for the growth rates A of the real infinitesimal perturbations tn,m and 
Tjn^rn [54] . Siucc the problem is Hamiltonian this matrix is infinitesimally sym- 
plectic and will have the simultaneous eigenvalues ±A, ±A* [55]. Due to the 
global phase invariance, A = is always a double eigenvalue. A stationary 
solution is linearly stable if all eigenvalues arc located on the imaginary axis, 
i.e., we can actually at most make a statement about marginal stability of the 
solution (meaning that a perturbation at most will grow with a rate that is 
polynomial in z). 

Extensive calculations, covering a large portion of parameter space, have been 
done for the most fundamental stationary solutions. In Fig. 1 the regions of 
stability (shaded) for the 1-site and the two 2-site solutions are shown for fixed 
norm J\f — 5. For the 1-site case there actually exist three different solutions 
in part of the region > 0.165 (the edge of the instability region), each ac- 
cessible by continuation but via different paths from the anti-continuous limit. 
Stability is indicated if one of the solutions is stable (there is no multistability 
of these three 1-site solutions and some details can also be found in Figs. 3 
and 4). The same behaviour appears also in ID [56], and is presumably an 
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(a) 2, (b) 2, (c) 2 




Fig. 2. Examples of solutions for Q2 = 0.2, Q3 = 0.5, Q4 = 0.25, Q5 = 0.05 
and J\f = 5. (a) 1-site solution with H = 19.7435 and A = 7.5019; (b) symmetric 
2-site solution with H = 19.9278 and A = 7.6357; (c) intermediate solution with 
n = 19.7310 and A = 7.5036. 

effect of the nonlinear coupling terms that can lead to a vanishing effective 
coupling between lattice sites. See further the discussion in Sec. 5. From Fig. 1 
we see that both stable and unstable solutions of all kinds exist and we note 
that the stability regions for the 2-site solutions together cover the region of 
instability of the 1-site solutions. For a quite large range of parameter values 
we thus have stable propagation of both types of solutions along the array of 
waveguides. The stability boundaries of the different solutions do not coincide. 
In the ID model [32,56], the stability boundaries very nearly coincide and in 
fact create a very narrow region of parameter space where the 1-site solution 
and one of the 2-site solutions have either simultaneous stability or instability. 
Moreover, the exchange of stability across this narrow range is connected to 
the existence of an intermediate asymmetric solution in this region, interpolat- 
ing between the two solutions. When crossing the region of stability inversion, 
the intermediate solution appears from a pitchfork bifurcation as one of the 
solutions (say the 1-site solution) changes its stability and then disappears in 
a reversed pitchfork bifurcation with the other (2-site) solution [32]. Here, in 
2D, although the stability boundaries lie rather far apart in parameter space, 
unstable asymmetric intermediate solutions still exist between the stability 
boundaries of the 1-site solutions and both 2-site solutions (Fig. Id) and ap- 
pear as the result of similar bifurcations. Examples of the solutions can be 
seen in Fig. 2. 

Moreover, comparing the Hamiltonians of the obtained solutions reveals the 
existence of zeros of the energy difference along some boundaries in param- 
eter space. These have been indicated in Fig. 1. To better clarify the rela- 
tions among the different solutions, bifurcation diagrams for Q4 = 0.25 and 
Q4 = 0.1 are shown in Figs. 3 and 4, with Hamiltonian Ti for the solutions 
plotted against the parameter Q^. Note that the solution that maximizes 7i 
always is stable. From Fig. 3 it is clear that the intermediate solution ap- 
pears in connection with the vanishing energy difference, i.e., the difference in 
Hamiltonian, of the 1-sitc and 2-site solutions ATi = Tii — 7^2, and through 
bifurcations provide the mechanism of their exchange of stability. The inter- 
mediate solution always lies very close but has a slightly lower value of 
and is therefore unstable, in contrast to the solutions in ID that may be sta- 
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-0.25 -0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 

Fig. 3. (Color online) Bifurcation diagram with Hamiltonian 7i versus the parameter 
Qs for Q2 = 0.2, Q3 = 0.5, Q4 = 0.25 and M = 5. 1-site solutions are indicated with 
circles (o), symmetric 2-site solutions with squares (□), anti-symmetric 2-site solu- 
tions with diamonds (O) and intermediate solutions with triangles (A). The markers 
are filled (unfilled) if the solution at that point is stable (unstable). Note also that a 
maximizer of the Hamiltonian always is stable. For some selected points (large cir- 
cles) a one-dimensional cross section along a lattice direction of the two-dimensional 
solutions is shown to give a view of their characteristic form. The intermediate solu- 
tions are always unstable and exist only in a limited parameter range. They emerge 
from pitchfork bifurcations and connect the 1-site solution with either of the 2-site 
solutions. Note that they appear in connection with the vanishing of the energy 
difference (see insets). 

ble [32]. From a mobility point of view, it is the energy of the intermediate 
solution that must be overcome to achieve mobility. This additional energy 
(the PN-barrier) is, with reference to the stable modes, small compared to the 
total energy (< 1%). Note also in Fig. 3 that for each 2-site solution there 
are two intersections with 1-site solutions and hence two zeros of the energy 
difference. One is, as expected from the energy minimum principle, present in 
the region of existence of the intermediate solutions, where the 1-site solution 
and the respective 2-site solution have simultaneous stability and are the sta- 
tionary modes with highest value of TC (see the zoomed areas). The energy 
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Fig. 4. (Color online) Same as in Fig. 3 but for Q4 = 0.1. Note here that 7i for the 
1-site solution is always larger than for the 2-site solutions. 

difference is also zero in a region where the involved solutions are not the 
ones with highest Ti, and in particular the 1-site solution is unstable in this 
region (see intersections at 7i ~ 15). This intersection, as can be seen from the 
solution insets in Fig. 3, is with a solution possessing the 'wrong' symmetry, 
i.e., the nearest neighbours of the centre site have the opposite sign compared 
to the symmetry or anti-symmetry of the respective 2-site solution. Therefore 
exchange of stability, or mobility, cannot be expected at this vanishing of the 
energy difference. 

In Fig. 3 we also see that the branches of the 1-site solutions terminate. Two 
of the branches will merge at ~ 0.01, thus the end of the branches are due 
to bifurcations. The behaviour of the three branches is not symmetric and the 
two branches lying close in the left part of the figure will not coincide. Instead, 
the upper branch of these two is terminated at ~ —0.14 in a bifurcation 
with a symmetric 5-site solution that at the anti-continuous limit has 5 in- 
phase neighbouring sites excited in a cross configuration. The lower branch 
can be continued to lower values of Qs where it will approach a plane wave 
solution. 



4 Mobility 



A bit surprisingly, when considering the results of the stability analysis of 
stationary solutions, the mobility of localized modes in the model (1) is very 
poor and in stark contrast to the very excellent mobility in the ID model [32] 
and the saturable model in ID [30] and 2D [31] near points of minimal PN- 
barrier. The method used to try to induce mobility is to apply a phase gradient 
^n,m ^ \E'„ me'*-*^'""'^*^'''"^ across one of the stationary real solutions. In appli- 
cations this can be achieved by launching the laser beam that excites the 
waveguides at an angle to the face of the array [5]. This perturbation will 
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Fig. 5. (Color online) Evolution of the centre of energy when perturbing a 1-site 
solution for Q2 = 0.2, = 0.5, (54 = 0.25, = 0.05 and = 5 (solution shown 
in Fig. 2a) with a phase gradient kx- In all cases W is sufficiently lowered to be below 
the value of the lowest stationary (the intermediate) solution with H. = 19.7310. The 
solid line is for = 0.1 {H = 19.7027), the dotted line for k^ = 0.2 {H = 19.5817) 
and the dash-dotted line for k^ = 0.5 {H = 18.7839). 

kick the solution in the direction of the phase gradient [24,30,31,32,36] and 
will have an effect similar to the marginal mode perturbation at an instability 
boundary described in [23,57], which in a complex formulation will correspond 
to a pure imaginary (velocity) perturbation to the real (position) solution. The 
phase gradient perturbation also has the nice feature of preserving the norm. 
Here we are strictly interested in axial propagation as this is the direction of 
possible propagation in the lattice suggested by the minimal PN-barrier of 
the stationary modes. Therefore ky = has consistently been used. Mobility 
in the diagonal direction would for the nearest-neighbour coupled lattice be 
connected to a vanishing energy difference for the 1-site solution and a 4-site 
solution (quadrupole). This does not appear in the present model and the 4- 
site solutions generally has a much lower value of Ti, than the 1-site solution. 
The reason that it is possible in the model in [31] is the saturable nonlinearity 
that will limit the amplitude of the excitations, bringing the singly excited 
and multiply excited solutions closer in energy for high enough power. The 
same effect is also the reason for the multiple zeros of the energy difference, 
ATC = 7^1 — 7^2, in saturable media, as more and more excited sites reach the 
saturation limit. For the model (1) we should expect at most one zero of the 
PN-barrier with fixed coupling parameters and varying norm A/". 

Trying to kick any of the stationary solutions in the region of existence of 
the intermediate solution, where the PN-barrier is minimal, will result in an 
oscillatory behaviour near the initial position. The centre of energy for the 
axial direction n, defined as 

'^CE = T7 H ^■^n,m , (9) 
■' ^ n,m 

will at most be displaced about one lattice site. The generic behavour is ex- 
amplified in Fig. 5 where the 1-site solution in Fig. 2a is perturbed with 
different phase gradients. Although the perturbation is large enough to over- 
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Fig. 6. Evolution of the centre of energy when perturbing an unstable symmetric 
2-site solution for Q2 = 0.2, Q3 = 0.5, = 0.25, = -0.06 and = 5 with a 
phase gradient kx = 0.2. The motion depends sensitively on initial conditions and 
boundary conditions. Here periodic boundary conditions on a 21 x 21 lattice is used. 

come the PN-barrier, set by the unstable intermediate solution, no mobility, 
except small displacements to neighbouring sites, is observed, indicating that 
the effect of pinning is stronger than suggested by the energy difference con- 
siderations. The final position of uqe is not even always in the direction of the 
phase gradient kick, and is located either on or between sites corresponding 
to one of the stable stationary solutions. Since the initial velocity of the ex- 
citation given by the perturbation can be large, the absence of mobility here 
is not the effect of resonance with the multi-spectral bands appearing at low 
velocity for travelling waves near a zero of the PN-barrier [53,58]. 

If we instead start with solutions away from the region of existence of the 
intermediate solutions only small oscillations around the stable stationary so- 
lution in the present region is observed. Displacement of the centre of energy 
a few sites is however possible. In Fig. 6 an unstable symmetric 2-site solution 
is perturbed and the centre of energy moves up to two sites. But the irregu- 
lar oscillations are due to a reconfiguration of the excitation from an initial 
more symmetric profile to, at the final position, a more anti-symmetric profile 
being closer to the stable stationary mode at the given values of the parame- 
ters. The motion is highly unpredictable and depends sensitively on the initial 
conditions and boundary conditions. 

At the edge of the stability boundary, where the bifurcation with the interme- 
diate solution occurs for the fundamental stationary modes, the eigenmodes 
related to the real eigenvalues that are responsible for the instability has a 
form that may be indicative of mobility. The real part of this growing mode 
is shown in Fig. 7 for the 1-site solution. Exactly at the stability boundary 
the imaginary part of the eigenmode is zero. Prom the form of the eigenmode 
we conclude that perturbing along this mode ('depinning mode' in [23]) will 
promote a change of shape of the 1-site solution towards a profile that resem- 
bles the symmetric 2-site solution (or the intermediate solution). The opposite 
relation holds at the stability boundary of the 2-site solution. Hence the form 
of the growing modes supports the heuristic view of mobility of narrow excita- 
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Fig. 7. Real part of the unstable eigenmode of the 1-site solution for Q2 = 0.2, 
= 0.5, (54 = 0.25, Qs = 0.0165 and J\f = b. By symmetry, the eigenmode is 
degenerate with a corresponding eigenmode oriented along the rh direction. The 
related eigenvalue is A = 0.4140. 
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Fig. 8. (Color online) Hamiltonian 7i, relative to the value of the 1-site solution (o), 
versus norm J\f for Q2 = 0.2, = 0.5, Qi = 0.25 and Q5 = 0.05. The intermediate 
solution (A) emerges from a pitchfork bifurcation with the symmetric 2-site solution 
(□) at ~ 2.1. It is always the solution with the slightly lower value of 7i and 
is unstable. As M increases, although very close in 7i, there will be no bifurcation 
with the 1-site solution, implying that the intermediate solution will not connect 
the 1-site and the 2-site solutions when considered as a function of the dynamical 
parameter M. The inset shows Ti versus N for the 1-site solution, as a reference. 

tions as a transformation between stationary solutions (see also Fig. 3 in [32]). 
But in the present model the stability boundaries are far apart in parameter 
space and such a transformation may not be possible. To be conclusive in 
this matter we will need to fix the coupling parameters and investigate the ex- 
change of stability as a function of the dynamical parameters Ti and M . If M is 
varied for a set of fixed coupling parameters inside the right region of existence 
of the intermediate solution in Fig. 3 we can find a situation when the inter- 
mediate solution will only bifurcate with the symmetric 2-site solution, and 
not with the 1-site solution for reasonable values of the norm {M < 100). This 
is illustrated in Fig. 8. In the event that a bifurcation occurs with both the 
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Fig. 9. (Color online) Hamiltonian 7i, relative to the value of the 1-site solution 
(o), versus norm Af for Q2 = 0.2, Q3 = 0.5, Q4 = 0.25 and Q5 = -0.14. The 
intermediate solution (A) appears at 1.4 and disappears at M ^ 10.9 in a 

bifurcation with the anti-symmetric 2-site solution (O). The inset shows 7i versus 
M for the 1-site solution, as a reference. 
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Fig. 10. (Color online) Hamiltonian TC, relative to the value of the 1-site solution 
(o), versus norm M for Q2 = 0.2, Q3 = 0.5, Q4 = 0.25 and = —0.1. In this 
case the intermediate solution (A) emerges from a pitchfork bifurcation with the 
1-site solution at ~ 2.4 and there is actually no zero of the energy difference. 
The anti-symmetric 2-site solution (O) always has a higher value of TC and is stable. 
The inset shows 7i versus M for the 1-site solution, as a reference. 

fundamental solutions involved (Fig. 9), the bifurcation points are far apart 
measured also in terms of the dynamical parameters. Measuring the distance 
between the bifurcation points in norm (AA/") we find that it is always compa- 
rable in size to the value of the norm at the bifurcation point with the higher 
value, i.e., AM /M ~ 1. The PN-barrier (A7i) is small compared to Ti for all 
values of the norm. Since the bifurcation points are either far apart or the in- 
termediate solution does not connect the 1-site and 2-site solutions at all, this 
indicates that there is no way to transform between the two. Put in another 
way, this suggests that there is no trajectory in phase space that passes close 
to the stationary configurations and consequently the mobility is very poor. 
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Our numerical observations lead us to the conclusion that mobility is resisted 
since the site-centred and bond-centred configurations are too far apart for 
the adjustment of M and Ti to be possible by exchange with an accompany- 
ing oscillating background, or in this case mainly an adjustment of H since 
the difference in energy is always relatively small. A tractable analysis could 
possibly be made along the lines of the method of an effective Hamiltonian 
where the dynamics is averaged over periodic trajectories in phase space [59]. 
The more explicit interpretation of the PN-barrier as referred to a trajectory 
in phase space would definitely be more suitable in the present context. How- 
ever, some extensions are necessary as an evaluation of the method shows that 
it gives an accurate description of mobility only close to limits where the trav- 
elling excitations are exact [59]. But this does not exclude that some insight 
might be given to the effect of pinning observed here. Further, as shown in 
Fig. 10, there are cases when the PN-barrier does not vanish at all, but the 
intermediate solution still exists. 



5 Compact solutions 

An interesting effect of the nonlinear couphng terms in Eq. (1) is that they 
may cancel the linear coupling and effectively give rise to a total coupling that 
is zero [32,42]. This property is also the criterion for the existence of exact com- 
pact solutions having strictly zero amplitude outside an interval [60,61]. The 
criterion will take the form of a parameter constraint and this implies that the 
compact solutions, as exact mathematical entities, are not robust with respect 
to parameter variation. They are interesting anyway from a physical perspec- 
tive, since a near perfect localization prevails in the neigbourhood of the exact 
parameter values, and this family of strongly localized solutions may very well 
be dynamically stable. The solution will develop an exponentially decaying 
tail, but the contrast, i.e., the ratio of the amplitude of neighbouring sites to 
the amplitude of the excited sites, is still large which can be desirable for appli- 
cations [24,32]. Further, although we speak of exact compact solutions in the 
context of the lattice equation (1), this does not imply compactness in the real 
system. The equation describes the amphtude of the modes in the respective 
waveguides, but the modes are themselves spatially exponentially decaying. 
The extreme localization represented by compact sohitions will nonetheless 
be preferable if several excitations propagate in the system at the same time 
as these will constitute excitations with least mutual interaction. The exis- 
tence of compact solutions has an impact on the properties of the system 
under study. 

The idea of the compactification is that for a particular amplitude the total 
coupling can vanish when the neighbouring sites have zero amplitude, leading 
to a decouphng of the lattice between these sites. Take, e.g., the 1-site exci- 



18 



tation with only the site (n, m) excited and all others at rest and put this 
ansatz into Eq. (8). From the equation on one of the neighbouring sites, e.g., 
site {n — l,m), we get Q2 + '2Q5\ipn,m\'^ = 0. For a compact solution with M 
excited sites, this constraint generalises to 

g2 + 2Q5^=0, (10) 

if all sites have the same amplitude. Note that < is required for the exis- 
tence of compact solutions (as Q2 > 0). This criterion corresponds to having a 
negative Kerr index of the nonlinear material surrounding the waveguides, or 
for BEC in an optical lattice an effective interatomic attraction in the dilute 
gas (as for ^Li). Solving for the 1-site and 2-sitc compact solutions the results 
are the same as in ID [32,42]. From Eq. (10) the constraint for the 1-site so- 
lution will correspond to the line Q5 = —0.02 in Fig. 1, that precisely cuts 
through the edge of the instability region of the 1-site solution and hence also 
seperates the two different regions of existence of the intermediate solutions. 
The loose interpretation of this constraint as an effective coupling is further 
strenghtened when considering that the two 2-site solutions will have the same 
value of Ti along the plane of existence of the 1-site solution. Note, e.g., that 
the lines of 7i intersect at Q5 = —0.02 in the bifurcation diagrams in Figs. 3 
and 4. 

Compact solutions with more than two sites excited can also be found. There 
are two different classes of real 3-site solutions, with the center site being 
either in-phase or out-of-phase with the two outermost sites. Having, in 2D, 
the freedom to choose the orientation of the three sites in the lattice as being all 
along one lattice direction or in an 'L'-configuration, the presence of coupling 
in two different lattice directions will also mean less freedom when compared to 
the ID system. The extra restrictions imposed for the 3-site solutions implies 
that all sites must have the same amplitude, which is not required in ID. 
Explicit calculations yield the same frequency A = 2Q2M /"^ for all solutions, 
and apart from the parameter constraint (10) also the constraint Q4 = ^Q^, 
where the upper sign is for the in-phase solution and the lower sign for the 
out-of-phase solution. A check of the stability shows that only the out-of-phase 
solution will be stable for J\f > 8.72, regardless of the orientation in the lattice. 
It is also possible to solve analytically for compact 4-site and 5-site solutions, 
some of which are stable for large enough values of the norm [M > 20). 

5.1 Vortices and complex solutions 

An interesting extension when going from one to two spatial dimensions is 
the possibhty to support localized excitations having a topological charge, 
which is proportional to the angular momentum carried by the excitation. 
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Fig. 11. (Color online) (a) The amplitude y^N'n,m = |^I'n,m| and (b) the phase 
Qn,m = arg{^'„_m} for a stable stationary compact vortex of topological charge 
S = 1 for M = 20, Q2 = 0.2, Q3 = 0.5, Qa = 0.04 and Q5 = -4Q2/AA = -0.04. 
The frequency is given by A = -3\/2Q2 + <53AA/4 + Q^M = 2.4515. 

We define the topological charge S as the number of complete 2tt twists of 
the phase when moving on a contour around the excitation. Such discrete 
vortex solitons have been observed in experiments [10,11] and described in 
the context of DNLS-type equations [62,63,64,65]. In the model (1) compact 
vortices can be supported. The constraint (10) must still be fulfilled, which 
means, as we have fixed Q2 and Qs by rescalings, that the parameter Q5 is 
fixed if the norm M and are used as free parameters. The most simple 
compact vortex is given by four excited sites in a square configuration with a 
phase difference of 7r/2 between neighbouring sites giving a topological charge 
S = 1. This solution will have the frequency A = Q^M /2 + Q^M and have 
a very small region of stability near (^4 = (e.g., for A/" = 10 it is stable 
for Q4 < 0.0023). For the square 4-site configuration placed diagonally in the 
lattice, with a zero-amplitude mediating site in the middle, the frequency is 
instead A = Q^M /2 and the region of stability is dramatically extended so 
that for A/" = 10 the configuration is stable for < 0.1634. The excited sites 
are not nearest neighbours in this case and will only weakly interact. In fact 
for A/" = 20 the instability will set in at Q4 = 0.165, which exactly corresponds 
to the instability of the compact 1-site solution for M = b. With more excited 
sites it is possible to find larger S = 1 vortices that are stable for certain 
parameter values. In Fig. 11 an example with eight sites is shown. Also higher 
order compact vortices exist. The configuration with eight sites on a square 
contour with phase difference 7r/2 between the sites will have charge S = 2 
and frequency A = Q^Af/A + Q/^M /2. It will also be stable for large norm 
(AA > 20) and small Q^^. (Note that, by contrast, this S = 2 configuration is 
always unstable in the standard on-site DNLS model [64].) 

The vortex solutions are necessarily complex solutions as they have a phase 
difference between neighbouring sites other than or tt, i.e., they cannot be 
transformed to real form using the global phase invariance of the DNLS-type 
equations. This implies that the norm current, Eq. (7), is non-zero and that 
there is a flow of norm around the contour of the vortex balanced according 
to the discrete continuity equation (6) with constant norm density N'n,m- As 
described for the ID model in [42] the presence of the nonlinear coupling 
will introduce a non-trivial phase dependence in the norm current density. 
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In particular, the norm current in Eq. (7) can become zero for specific 
values of the phase difference given by 

LiJbipn,m — i — ■ \^^) 

Thus under the constraint (11) stationary complex solutions can exist without 
having to form a closed loop for the norm current flow, which is the case 
for a vortex, simply because there is no current. The compact complex 2- 
site solution present in the ID model [42] exists also here, with cos (pn,m — 
Q2/2QJ\f and A = + Q^M - 3QI/2Q4A/". It is also stable for a range 

of phase twists, e.g., for A/" = 10 stability is found for 0.010 < Qa< 0.0420, 
equivalent to 0.2381 < cosipn^m ^ 1- Also compact complex solutions with 
more sites excited can be constructed, e.g., 4-site non- vortex solutions that 
generally seem to be unstable. Neither these solutions nor the vortex solutions 
above need to be compact, but can be continued to non-compact solutions 
by lifting the constraint set by Eq. (10). Although this may be an interesting 
further development the given examples serve to give an idea of the phenomena 
at hand. 



6 Power transport by plane waves 



Leaving the localized solutions, we will in this section focus on the transport 
properties of plane waves. The excitation of the lattice with a plane wave 
correponds to an even distribution of the energy over all sites in the lattice. 
A travelling wave on the form 

*n,m(^) = ^pe-'^^^+^+^^) (12) 

will have a norm density (power per site) Hn,m = P and be a solution of Eq. (1) 
for the frequency 

A= 2g3p + 2(Q2 + 8g5P)(cos(/? + cos0) 

+ 4Q4p(4 + cos 2^p + cos 2(f). ^ ' 

Note that the dispersion relation (13) can be derived from Eq. (5a) in the 
form A = dHn,m/dp. Though being a stationary solution, the plane wave (12) 
is associated with a balanced transfer of norm (power) between the sites due 
to the induced phase difference (wave number). This balance is described by 
the discrete continuity equation (6) with constant norm density. From the 
interpretation of ^7*^"^ = „ and J'^™-^ = J^"l^ as the norm transferred from 
a given site to its neighbours in the positive fi and m directions we can, 
from Eq. (7), obtain the norm current flowing in the lattice for a plane wave 
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solution, 



n,m 



dU 



dip 



-n 



n,TO ^ 

m 



2p sin V? (52 + 4(^4^ cos + 4(55P 
+ 2p sin (j) Q2 + '^Qip cos + AQ^p 



n 



(14) 



m. 



Further, it is convenient to introduce a polar representation of the current 
defined by the relation \J\e" = + ij("), so that 



J{n) 2 _^ Jim) 



(15) 



is the magnitude of the current and 



a — arctan 



' j-(m) ~ 



(±7r) e ] — TT, 7r[ 



(16) 



its direction measured counter clockwise relative to the n direction. Note that 
the power transported along the waveguides is proportional to p and that the 
quantities (14), (15) and (16) describe the transport across the array. For a 
finite array such transport can only be sustained if there is a supply of power 
at one boundary and a drain at the opposite boundary. 



6.1 Modulational stability 



To ensure that the plane wave (12) has stable propagation its behaviour under 
a modulation is investigated. Adding a modulating perturbation with wave 
numbers q and p we make the ansatz 



[Z C 



i(qn+pm) 



v*{z)e 



-i{qn+pm) 



X exp[—i{ifn + (pm + Az)] 



(17) 



in Eq. (1) [42,43,44,45,66]. Keeping only terms linear in u and v, the modula- 
tional stability can be determined from the equations 
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with 



2{Q2 + 8Q4P cos ip + SQ^p) sin ip sin q 

+ 2 ((52 + SQ^pcoscf) + 8Q5P) sin sin p , (19a) 

2QsP - 2Q2(cos ip + cos 4>) 
+ 4Q4p[2 cos q{l + cos 293) — cos 2ip] 
+ AQ4p[2 cosp(l + cos 2(f)) - cos 20] 

+ 2((52 + 8(55p)(cos(^ cosg + COS COS p) , (19b) 

2Q3P + 4:Q^p{cos 2(/9 + 2 cos g + cos 20 + 2 cosp) 
+ 8(55p[cos(^(l + cosg) + cos 0(1 + cosp)] . (19c) 



The plane wave wiU be stable if the eigenfrequencies u;± = a±V&^ — are real 
for all q,p — tt, tt]. By explicit calculation we find that the eigenfrequencies 
can be written on the form 



dp 



sm q + 



dp 



smp 
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4 sin^ f 



4sin2| 



X |/(<^)sin2! + /(0)sin2^ + p^ 



1/2 



(20) 



where we have introduced the (energetic) effective mass 



m. 



(n) 



1 a 
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dif^ p dip 

= -2((52 + 4(55P) COS if 



8Q4pcos2ip , 



(21) 



with a similar expression for m^\(f)), and the function 



/((^) = -2[{Q2 + I2Q5P) cos(^ + 4g4p(2cosV+l)], 



(22) 



entirely in correspondence with the result in ID [42]. From the first factor 
in Eq. (20) we see that if the effective masses have opposite signs there will 
be at least some values of q and p that will yield instability. An important 

observation that follows is that if or l/rriy^^ change sign due to tuning 

of some available parameters, the plane wave can only change its stability or 
remain unstable (unless the second factor in Eq. (20) changes sign at the same 
time) . 
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6.2 Power currents 



Generally, the direction of power flow in a lattice will be completely determined 
by the induced phase gradient of the plane wave. For the cubic on-site DNLS 
equation {Q4 — Q5 — 0), e.g., the direction of the norm current is fixed by the 
relative phases between the sites, tano; oc sin0/ sirup, while the magnitude 
is linear in the amplitude, oc p. Thus any control of the norm current 
by tunable parameters is restricted to the wave numbers of the plane wave. 
With current experimental techniques, a phase gradient is induced by having 
the exciting laser beam directed at an angle to the array of waveguides [5], 
so the wave numbers are easily adjustable but not rapidly tuned. Also, for a 
finite array the wave numbers will be quantized in order to satisfy the imposed 
boundary conditions. 

As discussed in [42], with the introduction of nonlinear coupling terms in 
the model (1), the norm current will have a nonlinear dependence on the 
amplitude. This will make p more suitable as a tunable parameter, also from 
an application point of view. Including the term, but keeping Q4 = (which 
would be relevant for apphcation to BEG in an optical potential [43,44,45]), 
will introduce a nonlinear dependence in the magnitude \^\ through the factor 
Q2 + '^QbP- The sign of this factor will determine in which of the two possible 
directions (a), difi'ering by an angle tt and fixed by the phase gradients if and 

0, the current will fiow. For < it is hence possible to make the current 
zero for a finite value of the amplitude p and at the same time reverse its 
direction, just as in ID [42]. However, as the factor Q2 + ^Qbp appears also 
in the effective mass (21), the change of sign will affect the stability of the 
travelling plane wave. In particular, as noted in Sec. 6.1, the plane wave cannot 
be modulationally stable for both directions of the current near the point of 
inversion. 

Inclusion also of the term will change this situation, as the zeros of j'^^^ 
and l/m^"* oc dj^"'^ /dip in general no longer will coincide. However, numerical 
investigations reveal that stable inversion is not possible in the axial direction, 

1. e., with = 0, TT, for any values of </? and Q4 > 0, but it is possible in other 
directions. Taking (/? = 7^ 0, tt will give a diagonal direction of the current, 
with a = 7r/4, —Stt/A depending on the relative signs of and J^"'\ and 
always stability near the point of inversion. An example of such a flip of 
direction as the amplitude of the plane wave is varied is shown in Fig. 12. An 
even greater flexibility is presented if different phase gradients are chosen for 
the n and m directions, (p ^ 4>. As illustrated in Fig. 12 we may have the 
remarkable situation that as the amplitude p is tuned the magnitude of the 
current is near constant while the direction is continuously changed an angle 
7t/2 over a stable range. Such an exact control of the power flow by simply 
changing the amplitude of the plane wave may flnd useful applications. One 
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Fig. 12. (Color online) (a) The magnitude 1^7"! and (b) the direction a of the norm 
current for varying amplitude p with the parameters Q2 = 0.2, Q3 = 0.5, Q4 = 0.1 
and = —0.1. The sohd line is for (/? = (/> = 7r/2, and the plane wave is stable for 
all p. The dashed line is for ip = 1.3 and 4> = 1-9) and the plane wave is stable for 
p > 0.119. 

possible application could be to use the array of waveguides, excited with a 
plane wave, as the link in a power coupling device. The advantage over a direct 
coupling would be a greater control of the flow facilitated by the control of 
the excitation in the waveguides. The effect of current inversion could form 
the basis for an amplitude-operated power-switch. However, any discussions in 
this direction other than mere speculative are beyond the scope of this paper. 



7 Conclusion 

In conclusion, we have studied the extension to two spatial dimensions of a 
one-dimensional model describing coupled optical waveguides going beyond 
the lowest-order approximation and including nonlinear coupling, which, e.g., 
is relevant for arrays of linear waveguides embedded in nonlinear media. The 
phenomenology of the ID model presented in a series of papers [32,42,56] has 
been revisited in the 2D setting. Especially, we have shown that a vanishing en- 
ergy difference between 1-site and 2-site stationary solutions and a subsequent 
small Peierls-Nabarro energy barrier, taking into account also the existence of 
asymmetric intermediate solutions, does not result in good mobility of highly 
localized modes. The main reason, and the difference to models that show 
good mobility [30,31,32], is that, despite small energy differences, the station- 
ary solutions are still in some sense far apart. Rather, the bifurcation points 
where the fundamental modes change their stability and where they exhibit a 
depinning mode promoting mobility are, although close in Hamiltonian (?-^), 
far apart in norm (A/") (Figs. 8 and 9). From this we conclude that for the sta- 
tionary solutions that have near equal values of both Hamiltonian and norm 
there is likely no trajectory in phase space passing close to them both. Thus, in 
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order to have good mobility of narrow modes in discrete systems, we need not 
only a small PN-barrier but also an exchange of stability between 1-site and 
2-site solutions taking place in the vicinity. We therefore also conjecture that 
the size of the oscillating background accompanying travelling modes [27,28] 
is not only related to the size of the PN-barrier but also to the difference in 
other conserved quantities (like Af) of the involved stationary modes measured 
at the points where their stability is changed. The lack of mobility near a zero 
of the energy difference between fundamental modes was observed also in [36], 
but no further analysis was carried out. 

Effects of a more delicate balance between linear and nonlinear coupling terms 
in Eq. (1) have also been discussed. These include the existence and stability 
of compact solutions, both discrete breathers and discrete vortex solitons. The 
latter have to our knowledge not been previously reported. Mathematically, 
these solutions are not robust with respect to parameter variation and require 
a balance of the coupling parameters in the equation for exact compactness, 
but near perfect localization persists in the neighborhood of the exact pa- 
rameter values. An interesting question is if they can prevail also beyond the 
tight- binding approximation of Eq. (1). The answer is likely that they are only 
exact solutions for the present model, but that they will correspond to modes 
of a higher degree of localization in more accurate models, i.e., they give an 
indication in which parameter regimes high localization can be achieved. As 
discussed in Sec. 5 the compact solutions do not represent a true compactifi- 
cation in the real system. 

Additional interesting effects present also in ID arise from a non-trivial depen- 
dence of the power flow in the lattice on the amplitude and phase difference of 
an excitation. Particularly, the current can become zero for a non-trivial phase 
twist which also in 2D leads to the existence of complex localized stationary 
modes with an open geometry, that may even be stable. Further, we have 
shown how the transversal flow of power in the array of waveguides can be 
controlled with great flexibility by excitation with plane waves, and in partic- 
ular explicitly demonstrated how the transport direction may be continuously 
tuned by amplitude variation. This may have applications for power-coupling 
devices. However, the use of a 2D array of the type studied here for multiport 
switching is discouraged due to the poor mobility of localized modes. For this 
purpose a ID array [32] or an array exhibiting a saturable [31] or quadratic [67] 
(see discussion below) nonlinearity instead of the Kerr nonlinearity is better 
suited. 

We end by commenting briefly on a number of works on related issues, that 
appeared in print after the original submission of our manuscript. 

In [68], AbduUaev et.al extended the Wannier-function approach of [18] for 
the continuous ID NLS equation with periodic linear potential to the case 
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with periodic variation also for the nonhnearity coefficient (which is the case 
also considered in our work). Within the tight-binding approximation, they 
derived under quite general conditions a lattice equation, which in a special 
case (cr = 1 in [68]) is equivalent to the ID version of our Eq. (1) as derived 
in [32]. An important conclusion of [68] is the crucial importance to include 
also the nonlinear coupling terms (corresponding to our Q4 and Q5), as they 
for specific choices of nonlinear interactions were shown to be comparable 
with, or even stronger than, the on-site nonhnearity. Thus, this supports the 
soundness of our approach to consider variation of the parameters Q4 and 
over a rather large range of values. 

Several works discussing properties of moving solitons and (vanishing) PN- 
barriers in ID generalized DNLS models have appeared. Dmitriev, Khare et 
al. [69,70] found analytically exact stationary and moving solutions for some of 
the exceptional, translationally invariant discretizations of [33,34,35]. Whether 
these are of relevance for any physically realizable system is, to our knowledge, 
still unclear. Pelinovsky et al. [71] developed a mathematical technique for 
analysis of persistence of traveling single-humped localized solutions from a 
certain limit, and found specifically that while travelling solutions terminated 
when continued from the integrable AL-limit of the Salerno model (corre- 
sponding to the development of a resonant tail [27]), they generally persisted 
in the translationally invariant models. Oxtoby and Barashenkov [72] used the 
method of asymptotics beyond all orders to evaluate the amplitude of radia- 
tion from a moving small-amplitude soliton in the saturable DNLS equation, 
and found it to be completely suppressed at certain 'sliding velocities', where, 
similarly as in [53], it was interpreted as an embedded soliton. The proper- 
ties of the travelling solitary waves in the saturable model were also analyzed 
numerically in more detail by Melvin et al. in [73]. 

For the two-dimensional case, Susanto et al. [67] studied the mobility of dis- 
crete solitons in a square lattice with quadratic (second-harmonic-generating) 
nonhnearity. In this case, due to the absence of collapse instability in the 
continuum limit, good mobility of stable solutions was observed as long as 
the coupling constants were not too small ('quasi continuum regime'). In this 
regime, mobility in arbitrary directions (not necessarily axial or diagonal) was 
observed, which is not surprising since discreteness effects are smoothened 
out for broad, continuum-like solitons. However, no non-trivial zeros of the 
FN barrier were found for this model, and consequently strongly localized so- 
lutions were reported to be immobile. Thus, in many aspects, the mobility 
properties of the 2D model with quadratic nonhnearity is similar to those of 
the ID cubic on-site DNLS model. 

Finally, a very recent preprint by Chong et al. [74] extended the analysis of 
the on-site cubic-quintic DNLS model of [38] to higher dimensions. As con- 
cerns the two-dimensional mobility, the results were shown to be very similar 
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to that of the saturable model [31]: enhanced mobility was found in regimes 
close to stability inversion and associated with appearance of asymmetric in- 
termediate solutions and low PN barrier. Thus, as far as we are aware, the 
present work still provides the only known explicit example of a model where 
these properties do not lead to an enhanced mobility. 
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